#
# Project: "Understanding Climate Change Conspiracy Theories 
#           And Their Believers: A Comparative Outlook"
#
# Authors: Jean-Nicolas Bordeleau & Daniel Stockemer (University of Ottawa)



######################### Load Packages and Dataset ############################

if(!require(tidyverse)) install.packages("tidyverse")
if(!require(dplyr)) install.packages("dplyr")
if(!require(plyr)) install.packages("plyr")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(ggeffects)) install.packages("ggeffects")
if(!require(sjPlot)) install.packages("sjPlot")
if(!require(modelsummary)) install.packages("modelsummary")
if(!require(QuantPsyc)) install.packages("QyantPsyc")


CCRS <- read.csv("CCRS Dataset (Final).csv")




################################# Analyses #####################################

### Table 1 ###

CTcrosstabs <- table(CCRS$ID_Q2, CCRS$CT2, useNA = "no")
prop.table(CTcrosstabs, 1)

FULLcrosstabs <- table(CCRS$CT2, useNA = "no")
prop.table(FULLcrosstabs)




### Figure 1 ###

formula_1 <- CT2 ~ ID_Q2 + ID_Q3 + ID_Q4 + ID_Q5 + ID_Q6 + ID_Q7 + ID_Q8 +
                    ID_Q10 + ID_Q20_7 + ID_Q23
model_1 <- lm(formula_1, data = CCRS, na.action=na.exclude)
modelsummary(model_1, stars = TRUE)

figure_1 <- plot_model(model_1, type = "std2", colors = "bw",
                       vline.color = "black",
                       rm.terms = c("ID_Q2 [Brazil]", "ID_Q2 [Canada]",
                        "ID_Q2 [Germany]", "ID_Q2 [Lebanon]",
                        "ID_Q2 [Morocco]", "ID_Q2 [South Africa]",
                        "ID_Q2 [United States]"),
                       sort.est=TRUE, title = "", 
                       axis.title = "Standardized Coefficients")

figure_1 + theme_bw() + ylim(-0.25, 0.25)




### Figure 2 ###

CCRS_AUS <- subset(CCRS, ID_Q1 == 4)
CCRS_BRA <- subset(CCRS, ID_Q1 == 8)
CCRS_CAN <- subset(CCRS, ID_Q1 == 1)
CCRS_DEU <- subset(CCRS, ID_Q1 == 3)
CCRS_ZAF <- subset(CCRS, ID_Q1 == 5)
CCRS_USA <- subset(CCRS, ID_Q1 == 2)
CCRS_LBN <- subset(CCRS, ID_Q1 == 7)
CCRS_MOR <- subset(CCRS, ID_Q1 == 6)

formula_2 <- CT2 ~ ID_Q3 + ID_Q4 + ID_Q5 + ID_Q6 + ID_Q7 + ID_Q8 + ID_Q10 + 
                   ID_Q20_7 + ID_Q23
model_AUS <- lm(formula_2, data = CCRS_AUS, na.action=na.exclude)
model_BRA <- lm(formula_2, data = CCRS_BRA, na.action=na.exclude)
model_CAN <- lm(formula_2, data = CCRS_CAN, na.action=na.exclude)
model_DEU <- lm(formula_2, data = CCRS_DEU, na.action=na.exclude)
model_LBN <- lm(formula_2, data = CCRS_LBN, na.action=na.exclude)
model_MOR <- lm(formula_2, data = CCRS_MOR, na.action=na.exclude)
model_ZAF <- lm(formula_2, data = CCRS_ZAF, na.action=na.exclude)
model_USA <- lm(formula_2, data = CCRS_USA, na.action=na.exclude)

all.models <- list()
all.models[[1]] <- model_AUS
all.models[[2]] <- model_BRA
all.models[[3]] <- model_CAN
all.models[[4]] <- model_DEU
all.models[[5]] <- model_LBN
all.models[[6]] <- model_MOR
all.models[[7]] <- model_ZAF
all.models[[8]] <- model_USA

figure_2 <- plot_models(all.models, m.labels = c("Australia",
                                                 "Brazil",
                                                 "Canada",
                                                 "Germany",
                                                 "Lebanon",
                                                 "Morocco",
                                                 "South Africa",
                                                 "United States"),
                        std.est = "std2",
                        colors = "Set1",
                        spacing = 0.6,
                        vline.color = "black",
                        legend.title = "Country",
                        axis.title = "Standardized Coefficients")
figure_2 + theme_bw() + ylim(-0.45, 0.40)



### Table 2 ###

EDUcrosstabs <- table(CCRS$ID_Q2, CCRS$ID_Q5, useNA = "no")
prop.table(EDUcrosstabs, 1)

GENDERcrosstabs <- table(CCRS$ID_Q2, CCRS$ID_Q4, useNA = "no")
prop.table(GENDERcrosstabs, 1)

CITYcrosstabs <- table(CCRS$ID_Q2, CCRS$ID_Q6, useNA = "no")
prop.table(CITYcrosstabs, 1)

IMMIGcrosstabs <- table(CCRS$ID_Q2, CCRS$ID_Q7, useNA = "no")
prop.table(IMMIGcrosstabs, 1)


### Table 3 ###

modelsummary(model_1, stars = TRUE)
lm.beta(model_1)


### Table 4 ###

modelsummary(all.models, stars = TRUE)


### Table 5 ###

variables <- c("CT2", "ID_Q3", "ID_Q4", "ID_Q5", "ID_Q6", "ID_Q7",
               "ID_Q8", "ID_Q10", "ID_Q20_7", "ID_Q23")

summary_table <- CCRS %>%
  select(all_of(variables)) %>%
  summarise_all(list(
    Mean = ~ mean(., na.rm = TRUE),
    SD = ~ sd(., na.rm = TRUE),
    Range = ~ diff(range(., na.rm = TRUE))
  )) %>%
  gather(key = "Statistic_Variable", value = "Value")

summary_table <- summary_table %>%
  separate(Statistic_Variable, into = c("Statistic", "Variable"), sep = "_") %>%
  spread(key = "Statistic", value = "Value")

kable(summary_table, caption = "Summary Statistics for Variables")

